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Abstract We employ an evolutionary algorithm to investi- 
gate the optimal design of composite protectors using one- 
dimensional granular chains composed of beads of various 
sizes, masses, and stiffnesses. We define a fitness function 
using the maximum force transmitted from the protector to a 
"wall" that represents the body to be protected and accord- 
ingly optimize the topology (arrangement), size, and mate- 
rial of the chain. We obtain optimally randomized granular 
protectors characterized by high-energy equipartition and the 
transformation of incident waves into interacting solitary pulses. 
We consistently observe that the pulses traveling to the wall 
combine to form an extended (long- wavelength), small- amplitude 
pulse. 
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1 Introduction 

One-dimensional (ID) lattices (chains) of particles interact- 
ing according to nonlinear potentials have been receiving in- 
creasing attention in the scientific community because of their 
special wave dynamics, which allows energy transport through 
solitary waves [39; 46; 59; 36; 9; 20; 21]. In the case of 
granular systems, particle interactions are strongly nonlin- 
ear because of nonlinear compressive forces and no-tension 
behavior [46; 59; 64; 28; 61]. As a result, granular lattices 
can support traveling solitary waves with compact support 
[46]. The evolution of nonlinear particle systems toward en- 
ergy equipartition (or thermalization), predicted by statisti- 
cal mechanics [ ], is also particularly interesting. Stable or 
transient energy transport through coherent modes (solitary 
waves, breathers, etc.) can develop ([20; 22; 23; 24; 18; 50; 
45; 17]), and eventual thermalization might not occur at all, 
as investigated in great detail for the Fermi-Pasta-Ulam prob- 
lem and related nonlinear lattice systems [9] . 



By designing protectors or containers optimally, the strongly 
nonlinear dynamics of granular systems can be exploited to 
produce fast decomposition of an external impulse into trains 
of solitary waves, energy trapping, shock disintegration, and 
more [29; 13; 47; 12;51; 15; 16; 64]. Furthermore, it has been 
emphasized that using a suitable randomization of the gran- 
ular system-involving, for example, the variation of particle 
sizes, masses, and materials-one might induce nonuniformity 
in the steady states in the velocity profiles; the appearance of 
negative velocities; marked thermalization; wave- amplitude 
decay; and anomalous features of wave propagation through 
interfaces between particles differing in masses, sizes, and/or 
mechanical properties [46; 13; 47; 29; 33]. 

When dealing with the optimal design of granular pro- 
tectors, one can optimize features such as particle distribu- 
tion, connectivity, size, and material through either discrete or 
continuous approaches. (These ideas are known, respectively 
as shape, topology, size, and material optimization.) Discrete 
approaches introduce suitable ground structure s-yN\\ic\\ refer 
to background structures in which the material densities of 
predefined connections are subject to optimization [54; 57; 
49; 7; 35]. Continuous models instead use homogenization 
theory, as they examine design domains with perforated com- 
posite microstructures [6; 5; 38; 31; 2; 7]. Well-established 
gradient-based optimization techniques include mathemati- 
cal programming, optimality criteria (that is, suitable math- 
ematical conditions defining an optimally designed structure) 
I ; ], sequential approximate optimization [63; 19; 37]. 
Available methods that are not based on gradients include 
simulated annealing [4; 62], biological growth [40; 42], and 
genetic and evolutionary algorithms [25; 32; 10; 53; 34]. 

Evolutionary Algorithms (EAs) provide a family of opti- 
mization methods inspired by Darwin's theory of evolution. 
They search for the best "phenotype" in a given population 
of candidate solutions by applying selection mechanisms and 
genetic operators similar to the intermingling of chromosomes 
in cell reproduction and replication. First, one evaluates each 
element (individual) of the population in terms of a quanti- 
tative fitness, which represents the feature that discriminates 
between phenotypes. One then mates the individuals using 



a recombination operator. Finally, one mutates a given per- 
centage of the individuals, thereby creating a new population. 
One then repeats these steps cyclically until some termination 
criterion is reached. The individual with the best fitness in 
the final population provides a guess of the global optimum 
of the fitness function. EAs are a natural fit for optimization 
problems in granular systems, as in such problems one can 
easily identify the system particles (i.e., the beads) with cells 
and their geometrical and mechanical properties (including 
radii, mass densities, elastic moduli, material types, and so 
on) with the corresponding genes. Furthermore, EAs require 
little knowledge of the search environment, can escape from 
local optima (in order to achieve a better global optimum), 
and are well- suited to problems with large and complex solu- 
tion spaces [ ] (such as those arising from the optimization 
of strongly nonlinear dynamical systems). 

The present work exploits EAs for the optimal design of 
composite granular protectors. We identify the fitness func- 
tion with the force Font transmitted from the protector to a 
"wall" that represents the body to be protected. We compute 
the performance of the candidate solutions under given im- 
pact loadings through a Runge-Kutta time-discretization of 
Hamilton's equations of motion. We adopt the hard- sphere 
model of interactions between adjacent beads for computa- 
tional reasons, as a very large number of simulations are re- 
quired by the optimization process. We also ignore dissipative 
effects, in accord with the standard models in the literature 
[46] . We note, however, that including relevant dissipative ef- 
fects [55; 56] such as friction, plasticity, large deformations, 
and so on, using (for example) time-stepping techniques of 
non-smooth contact dynamics [ ] or molecular dynamics 
[27], would not change the above EA framework. Dissipation 
is expected to enhance the effectiveness of the protector by 
further reducing the force amplitude transmitted at the wall, 
as shown (for example) in the experimental results reported 
below. 

In this paper, we investigate several optimization prob- 
lems. We focus, in particular, on topology, size, and material 
optimization of ID composite granular chains. We compare 
the dynamics of the optimized systems we obtain with those 
of granular protectors and special granular systems (sonic 
vacua) available in the literature [46; 16; 12; 29]. We show 
that the use of EAs offers a dramatic advantage in the design 
of granular protectors, leading to a significant decrease of the 
transmitted force. This EA-driven optimal design generates 
suitable topology, size, and material randomization by com- 
bining effects of wave disintegration and reflection at the in- 
terfaces between geometrical and/or mechanical discontinu- 
ities. A general feature we observe in the optimized protectors 
is the transformation of incident waves into a collection of 
interacting reflected and transmitted solitary pulses, which in 
particular form an extended (long- wavelength), small-amplitude 
wave that travels to the wall. We also find that optimization 
randomizes these systems (adding to their disorder) and pro- 
duces a marked thermalization. We constantly observe (in the 
absence of forced symmetry constraints) the appearance of 
soft/light beads near the wall, hard/heavy beads near the end 



impacted by the striker, and alternating hard and soft beads 
in the central section of the optimized chains. The observed 
"shock mitigation" behavior allows one to think of granu- 
lar protectors in a new way-as tunable kinetic systems rather 
than as purely dissipative systems. Consequently, they offer 
the exciting possibility of creating much more effective en- 
ergy transformation and shielding devices. 

The remainder of the paper is organized as follows. In 
Section 2, we formulate our mechanical and numerical mod- 
els. We then show how to optimize granular protectors in 
Section 3. We consider, in turn, topology optimization, size 
optimization, periodic sequences of optimized cells, and ma- 
terial optimization. As an extended example, we investigate 
the optimization of a container proposed by Hong [29] with 
both impulsive and shock- type loading. Finally, we summa- 
rize our results in Section 4. 

2 Mechanical and numerical modeling 

Consider a non-dissipative chain of N granular particles de- 
scribed by the Hamiltonian [46] 

^ /I 2 \ 
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where m^, qi and pi, respectively, denote the mass, the dis- 
placement from the "packed" configuration (particles touch- 
ing each other without deformation) and the momentum of 
the ith particle, Vi is the potential of the interaction force be- 
tween particles i and i + 1, and Wi is the potential of the 
external forces acting on the ith particle (including gravity, 
static precompression, etc.). We introduce an (TV + l)st par- 
ticle in order to model a wall that constrains the chain. In so 
doing, we assume that Pat+i = g'Ar+i = during the mo- 
tion. The Hamiltonian (2.1) yields a system of 2A^ first-order 
differential equations describing the motion of the system: 

dU dH 
Pi = Qi = ^ = 1,---,^, (2.2) 

oqi dpi 

to be solved with the initial conditions Pi{t = 0) = pf^\ 
qi(t = 0) = q^P , where t G [0,?] denotes the time vari- 
able, i indicates the final observation instant, and a dot over a 
variable denotes its derivative with respect to time. 

Assuming that stresses remain within the elastic threshold 
and that particle contact areas and velocities are sufficiently 
small, we introduce tensionless, Hertzian type power-law in- 
teraction potentials [46] 

= ^-a, [(g,-^,+i)+]^^+S (2.3) 
ni + 1 

where and are coefficients depending on material prop- 
erties and particle geometry, and (•)+ denotes the positive 
part of (•). Most of the examples that we examine in this pa- 
per are spherical grains, for which Hertz's law implies 

- 3f.+i(i-a + 3£:.(i-^f^,)' - 2 ' ^^-^^ 



where r^, 8i, and Vi denote, respectively, the radius, elastic 
(Young) modulus, and Poisson ratio of particle i. In the case 
of the granular container investigated by Hong [ ], we in- 
stead use the values shown in Table 1 for ai and n^. 
Additionally, let 
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and E = H = T , where T denotes the system's kinetic 
energy, V denotes the potential energy, and E denotes the 
total energy. We also introduce the local energies 
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Ei = + mqi-i - Qi) + Vi{qi - ^i+i)] (2.6) 

at each site (bead), and the energy correlation function (which 
is slightly different from that introduced in Ref. [17]) 



(2.7) 



where 
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(•) denotes the average over time (from to the current time), 
and 0) indicates how the energy is transferred between 
the different beads. Observe that 0) = corresponds to 
energy equipartition. 

Equations (2.2) can be solved numerically using a stan- 
dard fourth-order Runge-Kutta integration scheme (as dis- 
cussed in, for example, Ref. [h-o]) with a time integration step 
of 
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where Ci is the sound speed in the material for the ith parti- 
cle and k G (0, 1] is a scaling factor. Equation (2.9) gives a 
time-integration step of about 2 x 10~^ sfor/c = 0.1 and 1 
mm stainless steel bead chains (see Table 2), ensuring rela- 
tive errors lower than 10~^ in the total energy conservation 
for times up to few thousand /is. 

We now assume that the configuration of the granular sys- 
tem is described by a collection of design variables or "genes" 
(which can include particle radii, mass densities, elastic mod- 
uli, material types, etc.) 



subject to simple bounds of the form 



(2.10) 



(2.11) 



One can always assume assume that = and xf^ = 1, 
for alH G {1, M} through suitable rescaling of design 
variables. 




Fig. 2.1 (Color online) Diagrammatic representation of an evolutionary algorithm. 

Given an assigned x, a numerical simulation of the sys- 
tem dynamics under a prescribed impulse or shock loading 
gives the protection performance (fitness) / = ||Font||L°° 
of the corresponding design configuration. Here, Font de- 
notes the force transmitted from the system to the wall, and 
1 1 ^on^ II denotes its norm with respect to the Sobolev space 
L^([0,t]) [1]. The optimal design configuration Xopt can 
then be identified with the solution of the multivariate opti- 
mization problem, 

min f{x) , (2.12) 

which is expected to be influenced by multiple local optima. 
The problem (2.12) can be conveniently solved via EAs (see, 
for example, Refs. [25; 32; 26; 10; 53; 34]) through the cyclic 
iterative procedure illustrated in Fig. 2. 1. In the present paper, 
we will use the Breeder Genetic Algorithm (BGA) presented 
in Ref. [ ] . EGAs, in contrast to other EAs (in which the se- 
lection is stochastic), selects only from among the Tr% best 
elements of the current population of TV individuals (where 
Tr% denotes the so-called truncation rate) to be recombined 
and mutated (mimicking animal breeding). This feature makes 
the EGAs more efficient than standard EAs for performing 
optimization in large search spaces [44; 43]. 



3 Optimization of granular protectors 

We deal in the following sections with topology, size, and ma- 
terial optimization of ID composite granular protectors sub- 
ject to impulsive and shock- type loadings. We employ for- 
mula (2.9) with k = 0.1 for time discretization and always 
assume that genes are continuous variables ranging over [0,1] 
with a population size of 50 individuals; an initial, randomly- 
chosen, truncation rate (Tr) equal to 15%, Extended Inter- 
mediate Recombination (EIR) [ ], and a mutation rate in 
the interval [10%, 50%]. (We use the value 10% for size opti- 
mization, which consists of genuinely continuous genes and 
50% in all the other examples, which instead model discrete 
design variables using continuous genes.) EIR generates off- 
spring along the line defined by the parents in the search 
space and allows one to also create offspring outside of the 



segment joining the parents. See [14; 44; 43] for further tech- 
nical details of the employed BGA. 

The examples of Sections 3.1 and 3.2 consider chains 
of stainless steel beads, whereas those in Section 3.4 exam- 
ine a composite chain composed of polytetrafluoroethylene 
(PTFE) and stainless steel beads. We show the material prop- 
erties of these beads in Table 2. The final example, discussed 
in Section 3.5, considers a long composite chain-the pro- 
tector recently investigated by Hong [ ])-with the material 
properties shown in Table 1 . We studied protectors with one 
end in unilateral contact with a rigid wall (simulating the 
body to be protected) and the other end free. In most cases, 
we assumed that the free end was impacted by a striker; in the 
final example, we assumed that it was loaded by a prescribed 
force. We focused our attention on the short-term dynamics 
of the protector over an observation time slightly larger than 
that necessary to transmit the input actions to the wall. 
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matl 


2.0 


5657 


1.0 


mat2 


1.0 


5657 


2.0 


mat3 


0.3 


5657 


1.5 


mat4 


0.1 


5657 


1.5 



Table 1 Material properties (mass m, contact coefficient a, and con- 
tact exponent n) of the granular container investigated by Hong [29] 
(in abstract units, as discussed in the main text). 





p(kg/m^) 


8 (GPa) 


V 


stainless steel 


8000 


193.00 


0.30 


PTFE 


2200 


1.46 


0.46 



Table 2 Material properties (mass density p, elastic modulus and 
Poisson ratio v) of stainless steel and PTFE beads. 



3. 1 Topology optimization 

Nesterenko used the monicker sonic vacua to describe an un- 
precompressed (or weakly-precompressed) granular chains 
because the sound speed is zero or very small in such sys- 
tems [46]. He studied the behavior of two adjacent monodis- 
perse sonic vacua (2SV), characterized by a sharp variation 
in bead size (a stepped 2SV), under the impact of a striker. 
He observed two remarkable phenomena: disintegration of 
the incident pulse into a solitary wave train when it passes 
from the sub-chain with larger radius to the one with smaller 
radius; and a partial reflection in the opposite case (see also 
Ref. [33]). Here we examine the topology optimization of a 
stepped 2SV in order to determine the particle arrangement 
that minimizes Font under a given impact event. Figure 3.1 
shows different force-time histories in a 2SV hit by a striker 
at the sub-chain with larger radius. The system is composed 



of 20 large beads of radius r = vl = 3.95 mm, 20 small 
beads of radius r = rs = 2.375 mm, and the striker (par- 
ticle number 1), which has radius r = tl and initial veloc- 
ity v = 1 m/s (see § 1.6.10 of [ ]). The plots in Fig. 3.1 
show the force Fin at the contact between the striker and the 
first bead, the force F{i) that denotes the mean of the contact 
force between particles i and i — 1 and that between and 
i, and the force Font recorded at the wall. The observation 
time is 750 /is. All of the beads are made of stainless steel 
(see the material properties in Table 2). The F{i) plots for 
i > 21 in Fig. 3.1 clearly illustrate the aforementioned pulse 
disintegration phenomenon. One can also see that the fitness 
/= WFouth^ of the 2S Vis equal to 0.18 kN. 

We ran a topology optimization of the 2SV by introduc- 
ing M = TV = 40 genes Xi related to the radius size (large 
or small) of the different beads. (This does not include the 
striker-particle number 1 -which is assumed to have a large 
radius.) We defined the genes so that Xi G [0, 0.5] implies 
r^+i = rs, whereas Xi G (0.5,1] implies r^+i = r^. We 
used a penalty technique to constrain the number of parti- 
cles with large and small radii to each be equal to 20; that 
is, we assigned a very large fitness / to (unfeasible) solutions 
that do not satisfy this criterion. We show the BGA-optimized 
system and the corresponding force-time plots in Fig. 3.2. 
The optimized system has many large beads near the end of 
the chain that is hit by the striker (shown on the right), small 
beads near the wall, and an alternation of sequences of mul- 
tiple consecutive large and small particles in the center of the 
system. We obtained a stable solution (i.e., a solution with 
constant best fitness) after about 340 generations of the al- 
gorithm. Observe that pulse disintegration appears very early 
(within the first few beads) in the optimized system, so that 
the leading solitary wave transforms into a train of interact- 
ing, small-amplitude pulses. This configuration exhibits a fit- 
ness of about 0.049 kN, which is almost four times smaller 
than that of the 2SV. 

We compare the energies (as a function of time) of the 
stepped 2SV and the optimized system in Fig. 3.3 over a time 
window preceding the achievement of a loose state (in which 
there are no interaction forces).^ We obtained this by restrict- 
ing the energy time-histories up to the first instant t > for 
which T = 0.99£;. The kinetic energy T of the 2SV shows 
a marked peak when the leading wave passes from the larger 
sub-chain to the smaller one and valleys when the wave is 
reflected at the wall. The potential energy V behaves in the 
opposite manner because the total energy is conserved. In the 
optimized system, on the other hand, the valleys and peaks 
of T and V arise earlier during wave propagation, and the 
peaks of the potential energy remain markedly lower than 
those observed in the 2SV. Denoting by (T) and (V) the 
time-averaged values of T and V, respectively, over a time 
of 1000 ms from the striker impact, we observe that in both 
the 2SV (which has {T)/{V) ^ 1.55) and the optimized sys- 

^ In the absence of gravity and precompression, a loose state is 
reached after a sufficiently long time because the granular chain is 
constrained only at one end. The dynamics evolve so that the inter- 
actions go to zero. 
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Fig. 3.1 (Color online) Force versus time plots for the stepped two sonic vacua (2SV). 
The striker impacts the end with larger-radius beads. 
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Fig. 3.2 (Color online) Force versus time plots in the topology-optimized system. 
(Compare this to the (unoptimized) stepped 2SV configuration in Fig. 3.1.) 



tern (which has {T)/{V) ^ 1.86), the ratio {T)/{V) deviates 
from the value 1.25 predicted by the virial theorem of sta- 
tistical physics [ ]. Nesterenko observed similar results us- 
ing randomized granular chains subject to piston-like impacts 
[46]. 

Figure 3.4 shows density plots of particle energies Ei for 
the stepped 2SV and the optimized system. In each plot, the 
horizontal axis indicates the particle site, the vertical axis 
shows the progressive step number (we produced a plot for 
every five integration steps), and the shading gives a den- 
sity plot of the energy Ei normalized to unity (i.e., the en- 
ergy divided by its maximum value, among all the beads). 
One can clearly recognize the impulse disintegration phe- 
nomenon in the stepped 2SV when the incident pulse passes 
from the large-bead regime to the small-bead regime. In this 
system, pulse reflection occurs only at the wall (and not along 
the body of the chain) during the first transmission. Note 
that when already-reflected pulses pass from the small-bead 
regime to the large-bead regime, they are reflected for a sec- 
ond time. In the optimized system, one observes a combina- 
tion of disintegration and reflection of traveling pulses along 
the entire chain. One also observes the production of interact- 
ing pulses that travel in opposite directions. In particular, the 
pulses moving toward the wall combine to form an extended 
(long- wavelength), small- amplitude wave that is clearly visi- 
ble in Fig. 3.2. Figure 3.5 shows the time histories of the en- 
ergy correlation function for the stepped 2SV and topology- 
optimized systems, revealing that the latter exhibits a faster 
and stronger thermalization (i.e., equipartition of energy) than 
the former. 



3.2 Size optimization 

Doney and Sen recently studied the energy absorption ca- 
pabilities of ID granular protectors consisting of "tapered" 
and/or "decorated" chains [15; 16]. The simplest type of ta- 
pered chain is composed of a sequence of progressively larger 
or progressively smaller beads, and a decorated chain is a 
composite chain obtained by placing interstitial small grains 
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Fig. 3.3 (Color online) Energy versus time plots in the stepped 2SV and in the 
topology-optimized system (energies are in mJ, and times are in /x s). 
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Fig. 3.4 (Color online) Density plots of particle energies normalized to unity. Horizon- 
tal axes are labeled according to particle site and vertical axes give the time step. 
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Fig. 3.5 (Color online) Energy correlation function versus time for the topology- 
optimized and stepped 28 V systems. 



between the large grains in a monodisperse chain. Using an- 
alytical and numerical techniques, Doney and Sen observed 
marked energy absorption in highly tapered chains. For their 
analytical work, they employed the hard- sphere approxima- 
tion. In their numerical studies, they utilized hydrocode sim- 
ulations up to very high impact velocities (up to 1 km/s). In 
Figs. 3.6 and 3.7, we show several numerical force recordings 
that correspond, respectively, to a decorated and a tapered 
chain impacted by strikers with different velocities. In both 
cases, the initial force peak is Fin = 1.4 kN. The decorated 
chain consists of a sequence of dimers formed by alternating 
r = tl = 3.243 mm and r = rs = 0.973 mm stainless steel 
beads with total length 7.78 cm (not counting the striker). 
The striker (3.243 mm) impacts this chain with an initial ve- 
locity of 18.88 m/s. The tapered chain, on the other hand, is 
composed of 20 stainless steel beads with decreasing radii, 
ranging from r = tl = 5 mm (the striker, which impacts 
the chain at a speed of 7.20 m/s) to r = rs = 0.675 mm 
over a length of 7.78 cm (not counting the striker). That is, 
the tapering ratio is Qs = Ti+i/r^ =0.1. The fitness of the 
decorated chain is about 1.7 kN, whereas that of the tapered 
chain is equal to about 0.75 kN. 

We carried out a size optimization of the above systems, 
introducing M = 19 genes Xi related to the radii of the dif- 
ferent beads (not including the striker), which in the present 
example (in contrast to the previous one) were assumed to 
change continuously in the interval [rs^VL] (with r^+i = 
Xi{rL — rs))- We ran a BGA optimization with a striker 
radius always equal to 5 mm and an impact speed equal to 10 
m/s. We also constrained the total length of the system to re- 
main equal to 7.88 cm (not including the striker). We obtained 
constant best fitness and the solution shown in Fig. 3.8 after 
about 590 generations. The fitness of the size-optimized sys- 
tem (0.42 kN) is about 1.8 times smaller than that of the (un- 
optimized) tapered chain. Observe that there is simple reflec- 
tion at the wall in the decorated chain (see Fig. 3.6), signif- 
icant pulse disintegration in the tapered chain (see Fig. 3.7), 
and a transformation of the leading solitary pulses into an ex- 
tended (long- wavelength), small- amplitude wave in the op- 
timized chain (Fig. 3.8). The choice of the fitness parame- 
ters used in this study differs from Doney and Sen's [ ], 
which instead minimizes the kinetic energy ratio Kout / Kin 
between output and input. Because of the continuum formu- 
lation of the genes in the current examples, our work encom- 
passes all geometries considered in Ref. [ ' We compare 
the energy profiles of the decorated, tapered, and optimized 
chains in Fig. 3.9. We observed the highest (T) / (V) in the ta- 
pered chain ((T) / (V) ^ 1.94), resulting from the anticipated 
evolution of this system toward a loose state. Figure 3.10 
shows the density plots of particle energies for the three sys- 
tems under examination. As in the previous example, one can 
clearly observe from the plots the effects of combined wave 
disintegration and reflection in the optimized system. The 
profiles of the energy correlation function shown in Fig. 3.11 
indicate that the size-optimized and the tapered chains both 
evolve toward thermalization, with slightly faster decay of 
C(t,0) in the former system. We show in the next section 
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Fig. 3.6 (Color online) Force versus time plots in a decorated chain. 
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Fig. 3.7 (Color online) Force versus time plots in a tapered chain. 
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Fig. 3.8 (Color online) Force versus time plots in the size-optimized system. 

that a similar behavior can also be induced in a long dimeric 
system (i.e., in a long decorated chain) by introducing suit- 
able alterations of the periodic particle arrangement (i.e., by 
introducing another form of disorder into the system ). 



3.3 Periodic sequences of optimized cells 

We now consider periodic sequences of the 19-particle size- 
optimized cells in the decorated chain shown in Fig. 3.8. As 
discussed above, the single optimized cells that we obtained 
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Fig. 3.9 (Color online) Energy versus time plots in the decorated, tapered, and size- 
optimized chains (energies in J, times in /i s). 
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Fig. 3.10 (Color online) Density plots of particle energies normalized to unity. Hori- 
zontal axes are labeled according to particle site and vertical axes give the time step. 
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Fig. 3.11 (Color online) Energy correlation function versus time for the examined sys- 
tems. 



can be viewed as disordered configurations, so it is interest- 
ing to investigate the effects on the wave dynamics of periodi- 
cally repeating such a structure to obtain a "quasi-disordered" 
configuration. As can be seen in Fig. 3.12, a reasonably local- 
ized wave structure does develop as long as there are enough 
periods, just as with periodic arrangements of simpler cells 
[51; 52]. However, the original wave disintegrates and emits 
a significant number of secondary pulses, so that a stable co- 
herent structure does not form. 

Using long-wavelength asymptotics, one can obtain a non- 
linear partial differential equation (PDE) description of the 
decorated chain in the continuum limit [46; 5 1 ; 52]. This PDE 
has known compact solitary wave solutions, which we illus- 
trate in the top panels of Fig. 3.13. However, adding even a 
small number of impurities to the system can change things 
completely (introducing some pulse disintegration, pulse re- 
flection, and thermalization), although one still obtains a basi- 
cally localized pulse. The impurities we consider here consist 
of particles of radius 5 mm and mass 2.31 g, so that they are 
much larger and heavier than the other beads in the chain. 
Throughout the region of the chain that has impurities, we 
place one of them every 19 particles, giving a cell length that 
is the same as that in the periodically repeated size-optimized 
chain of Fig. 3.8. The second through fourth rows of Fig. 3.13 
contain regions of different lengths that contain these peri- 
odic impurities. In each case, the last impurity is placed be- 
fore bead 1000. The first impurity is in particle 501 in the 
second row of Fig. 3.13 The third row of the figure is for a 
chain with impurities every 19 beads starting from bead 801, 
and the bottom row is for a chain with impurities every 19 
beads starting from bead 95 1 (so that there are three impuri- 
ties in total-at beads 970 and 989-in this last example). As 
shown in these plots, the insertion of such heavy impurities 
leads to partial reflections of the wave, some thermalization, 
significant pulse disintegration, and even to a bit of trapping 
(see the bottom left panel). Also observe in the bottom row 
that we have induced delays in the wave reflection. By tun- 
ing the material properties carefully, one can perhaps opti- 
mize the properties of such wave trapping so that they can 
be exploited in applications. Moreover, these numerical ex- 
periments also illustrate the much more complicated series 
of secondary pulse emission and partial wave reflections that 
occur in the periodic sequence of optimized/randomized cells 
inRef. 3.12. 
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Fig. 3.12 (Left) Density plot (colored by force, with larger values given by darker shad- 
ing) of the optimized chain in Fig. 3.8. (Right) Force (in kN) versus time plot for particle 
1500 of this chain. 
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Fig. 3.13 (Left) Density plots (colored by force, with larger values given by darker 
shading) for (top row) decorated chain of Fig. 3.6 and for decorated chains with a single 
impurity in each 19-particle cell between beads 501 and 1000 (second row), 801 and 
1000 (third row), and 951 and 1000 (bottom row). (Right) Force (in kN) versus time 
plots for particle 1500 in each of these chains. 



The study of soliton-like pulses in perturbed uniform Hertzian 
chains was reported earher [60; 30], suggesting their use as 
possible systems to detect buried impurities via the analy- 
sis of back- scattered signals. Our results show that similar 
phenomena can also be observed in the "quasi-disordered" 
systems; as shown in Fig. 3.13, the mass and position of the 
defects in the chain detectably shift the reflection and the ra- 
diation. 

It would be interesting to extend this type of discussion 
by considering increasingly disordered configurations, such 
as systems with quasiperiodic arrangements of cells (with 
various lengths and component particles) rather than periodic 
ones. Some preliminary research in this direction (using, for 
example, arrangements that follow Fibonacci sequences) has 
appeared recently in the literature in order to study Anderson 
localization in atomic chains [ ; ] . It would similarly be 
interesting to consider systems with random arrangements of 
cells. 



3.4 Material optimization 

In Ref. V^], Daraio, et al. investigated the optimization of a 
composite granular protector consisting of 22 stainless steel 
beads and 10 PTFE beads (see Table 2) with a uniform radius 
of 2.38 mm. The authors examined different design solutions, 
based on material distribution, using both numerical and lab- 
oratory experiments. The protector they considered was im- 
pacted by an AI2O3 striker (0.47 g) with initial velocity of 
0.44 m/s and was initially precompressed by a static force 
of 2.38 N. Using piezosensors embedded in selected parti- 
cles, they obtained laboratory measurements of force versus 
time profiles in several beads and compared them against nu- 
merical predictions. Here we report analogous experiments to 
confirm the results obtained from the optimization analysis. 
We used a four-garolite-rod stand as the holder for the beads 
and assembled sensors as described in Refs. [51; 52]. We se- 
lected a steel particle (0.45 g) as the striker and recorded the 
traveling signal with a TKTDS 2024 oscilloscope (Tektronix, 
Inc.). The sensors (PiezoSystems, Inc.) were calibrated by 
conservation of linear momentum. The 1.86 N static precom- 
pression included the preloading of the topmost particle with 
about 190 g of symmetrically suspended masses. 

Figure 3.14 shows some numerical and experimental force 
recordings in a soft-hard- soft configuration (see Fig. 2 of Ref. [12]) 
with two sequences of five PTFE beads at both the top and 
bottom of the chain. This system had the minimum value of 
Font of all of the configurations considered in [ ]. (A differ- 
ent hard- soft-hard- soft-hard configuration minimized Font / 
Fin ratio, as shown in Fig. 3 of Ref. [12].) Observe the good 
qualitative agreement between numerical and experimental 
results over the initial phase of the pulse propagation. The dy- 
namics of the experiments and numerics subsequently deviate 
from each other, as the laboratory tests are affected by dissi- 
pative effects (not included in the numerics). The latter can 
arise from effects such as friction, inelastic collisions, vis- 
cous drag, etc. The experimental traveling wave is thus pro- 
gressively damped as it travels through the chain, resulting in 
an even better protector. 

We carried out a material optimization by introducing 32 
genes Xi related to the material identification of the individ- 
ual beads (not including the striker). They are defined such 
that Xi e [0,0.5] implies that the {i + l)th bead is made of 
PTFE, whereas Xi G (0.5,1] implies that the same bead is in- 
stead made of stainless steel. (We assumed that the striker was 
always made of steel.) We constrained the total number of 
PTFE beads to be equal to 10 through a penalty technique. We 
also introduced an additional gene (so that the total number of 
genes M is equal to 33) related to the intensity of the preload, 
allowing the static precompression force Fq to vary continu- 
ously within the interval [0, 2.38] N (Fq = 2.88x33 N). The 
optimized system, obtained after about 110 BGA generations, 
is shown in Fig. 3.15 together with the corresponding numer- 
ical and experimental force plots. As in the previous exam- 
ples, observe that the material-optimized system contains soft 
beads near the wall, hard beads near the end impacted by the 
striker, and alternating hard and soft beads in the central sec- 
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Fig. 3.14 (Color online) Force versus time plots in a soft-hard-soft granular chain. (The 
applied precompression was added to the force profiles experimentally recorded though 
piezosensors.) 
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Fig. 3.15 (Color online) Force versus time plots in the material-optimized system. (The 
applied precompression was added to the force profiles experimentally recorded though 
piezosensors.) 



tion of the chain. We computed the optimal precompression 
force to be about 1.86 N. Both the numerical and the experi- 
mental force plots of Fig. 3.15 show that the leading solitary 
wave first decomposes into a train of small pulses and subse- 
quently mutates into an extended (long- wavelength), small- 
amplitude wave. The density plots of Fig. 3.16 illustrate the 
mechanisms of wave disintegration and reflection character- 
izing the dynamics of these systems. The experimental results 
presented in this article and in Ref. [12] confirm the enhanced 
performance of the BGA-optimized system (minimum Fout), 
as compared to all the other examined protectors. 



Soft- hard-soft svstem 



Optimized system 




5 10 15 20 25 30 



Fig. 3.16 (Color online) Density plots of particle energies normalized to unity. Hori- 
zontal axes are labeled according to particle site and vertical axes give the time step. 



3.5 Long composite protector 

In a recent paper [ ], Hong investigated a long ID compos- 
ite granular protector (or "energy container") consisting of 
nine 20-bead sections. The beads in this chain were made of 
four different materials with varying particle mass m, contact 
stiffness a, and contact exponent n (see the parameter values 
in Table 1). Hong used abstract units, introducing factors of 
10-^ m, 2.36 X 10-^ kg, and 1.0102 x 10"^ s to convert 
the adopted units of length, mass, and time, respectively, into 
real units. The terminal and central sections of the protector 
are composed of a (linear) material ("material 1") character- 
ized by a contact exponent n = 1 and mass m = 2. The inner 
sections are composed of beads (made of nonlinear materials) 
that have different masses and contact exponents greater than 
1 (materials 2, 3, and 4). This is used to simulated sharp con- 
tact surfaces and rough materials such as sand (see Fig. 3.17). 
The container has a reflection symmetry about its center and 
the initial distance between particle centers of mass is uni- 
formly equal to 200 along its body. 

Hong analyzed the behavior of this container (and vari- 
ations thereof) by considering its dynamics after the impact 
of a striker of mass m = 100 traveling with speed v = 10. 
He ran numerical simulations of the chain dynamics, employ- 
ing a hard- sphere model and introducing lateral constraints 
through additional beads consisting of very heavy grains (m = 
100). He found a universal power-law scaling for how long 
it took the energy to leak from the protector to the lateral 
sides. That is, the energy remaining in the protector is given 
by Er = At~^, where A is a constant that depends on the 
protector construction, t is the time, and 7 is a constant (that 
Hong estimated in Ref. [29] to be about 0.7055) independent 
of the protector construction. 

We carried out a joint topology-material optimization of 
Hong's container (shown in the top panel of Fig. 3.17), in- 
troducing M = 90 genes Xi that characterize the material 
identification of each individual bead. The conditions Xi G 
[0,0.25],Xi e (0.25, 0.5], G (0.5, 0.75], andx^ G (0.75,1] 
respectively imply that the ith bead of one half of the protec- 
tor is composed of material 1, 2, 3, and 4. We enforced the 
symmetry with respect to the center of the chain by suitably 
relating the material identification numbers of the 90 grains of 
the second half to those of the grains in the first half. We did 
not enforce any constraints on the numbers of beads of the 
different materials. Due to the central symmetry constraint, 
the optimized protector is not allowed to have different con- 
structions near the impacted and the constrained ends, in con- 
trast to the protectors we examined in the previous sections. 

3. 5. 1 Impulsive loading In our first optimization procedure, 
we considered the impact of an m = 2 striker (material 1) 
traveling with speed v = 10. We show the corresponding op- 
timized protector, which we obtained after about 400 EGA 
generations, in Fig. 3.17. This optimal impulse absorber has 
nonlinear beads near its extremities and in its center and se- 
quences of linear beads in its remaining sections (a few of 
them are also near the center of the half chain). We show the 
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Fig. 3.19 (Color online) Energy versus time plots in Hong's container and the opti- 
mized chain under impulsive loading. 



Optimal shock absorber x 1/2 
matlae 

mat2io 
matS^i 

mat42i 



Fig. 3.17 (Color online) Hong's container and optimized composite long chains (show- 
ing half of the chain; the other half is obtained by reflection symmetry). 
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Fig. 3.18 (Color online) Force versus time plots in a long composite chain subject to 
impulsive loading (force values divided by 1000). 



Fig. 3.20 Density plots of particle energies normalized to unity. Horizontal axes are 
labeled according to particle site and vertical axes give the time step. 
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Fig. 3.21 (Color online) Energy correlation function versus time for the examined sys- 
tems. 



corresponding force plots and energy profiles (as well as the 
ones for Hong's container) in Fig. 3.18 and 3.19, respectively. 
Observe that the optimized scheme transmits to the wall a 
maximum force (126) that is about three times smaller than 
that transmitted by the basic scheme (371). As in the previ- 
ous examples, the initial pulse is progressively disintegrated, 
reflected, and transformed into an extended wave within the 
optimized system (as confirmed also by the density plots of 
Fig. 3.20). Hong's container instead shows wave reflection 
only when the incident wave crosses the central section of the 
system. The time histories of the energy correlation function, 
depicted in Fig. 3.21, highlight the fact that the optimized 
system spreads out energy (i.e., thermalizes) on a faster time 
scale than Hong's container. 

3.5.2 Shock-type loading We also carried out an optimiza- 
tion procedure using a shock-type force profile with constant 



intensity F = 1000 on the first bead (composed of mate- 
rial 1), external to the container, for a time equal to 0.25. 
We show the corresponding optimal container, which we ob- 
tained after about 200 BGA generations, in Fig. 3.17. The 
force-time plots of this optimal shock absorber and those 
of the basic Hong scheme, with the shock-type loading, are 
shown in Fig. 3.22. Observe that the basic protector transmits 
to the wall a peak force (1300 units) larger than the applied 
shock, whereas the optimized system is able to reduce the 
shock at the wall up to 770 units (a 23% reduction). Note 
additionally that the optimal shock absorber is characterized 
by heavy, linear (material 1) grains near the extremities (see 
Fig. 3.17). This is likely due to the symmetry constraint dis- 
cussed above. Figure 3.22 also shows that the input shock 
gets weakened when traveling along the optimal protector. 
We show the density plots of particle energies in these two 
systems in Fig. 3.23. 
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Fig. 3.22 (Color online) Force versus time plots in a long composite chain subject to 
shock-type loading (forces values divided by 1000). 
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of impulsive loading and one that had heavy, linear beads to- 
ward the ends in the case of shock- type loading. 

The present research paves the way for many interesting 
developments, as our approach can be generalized to numer- 
ous situations. First, the techniques we employed can be ap- 
plied to more intricate experimental configurations-including 
two-dimensional systems, three-dimensional systems, systems 
composed of ensembles of particles with non-spherical ge- 
ometries or even layered materials. Second, one can incorpo- 
rate additional physical effects, such as dissipation and more 
complicated contact mechanics. Third, one can generalize the 
methods themselves by, for example, adopting continuous op- 
timization techniques such as the material distribution method 
[7] and the formulation of multiple- scale approaches that in- 
volve scale-dependent interaction forces. 
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Fig. 3.23 Density plots of particle energies normalized to unity. Horizontal axes are 
labeled according to particle site and vertical axes give the time step. 



4 Conclusions 



In summary, we used an evolutionary algorithm to investi- 
gate the optimal design of composite granular protectors us- 
ing one-dimensional chains of beads composed of materials 
of various size, mass, and stiffness. Identifying the maximum 
force Fout transmitted from the protector to a "wall" that rep- 
resents the body to be protected as a fitness function, we op- 
timized the topology (arrangement), size, and material of the 
beads in the chain in order to minimize Fout- We considered 
several examples that were investigated recently in the lit- 
erature, including stepped two sonic vacua, tapered chains, 
decorated chains, and a recent configuration due to Hong. 

The optimization procedure, driven by a Breeder Genetic 
Algorithm, produced (optimally) randomized/disordered sys- 
tems, along which the incident waves were disintegrated and 
reflected, exhibiting marked thermalization. Additionally, the 
solitary pulses traveling to the wall combined to form ex- 
tended (long- wavelength), small- amplitude waves. In the ab- 
sence of enforced central symmetry, the optimal configura- 
tions had soft/light beads near the wall, hard/heavy beads 
near the loaded end, and alternating hard/heavy and soft/light 
beads in the remaining part of the chain. In the presence of 
central symmetry, we instead obtained an optimal configura- 
tion that had light, nonlinear beads toward the ends in the case 
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